// ****************************************************************************** 
// SEEDS Superpixels
// ******************************************************************************
// Author: Michael Van den Bergh
// Contact: vamichae@vision.ee.ethz.ch
//
// This code implements the superpixel method described in:
// M. Van den Bergh, X. Boix, G. Roig, B. de Capitani and L. Van Gool, 
// "SEEDS: Superpixels Extracted via Energy-Driven Sampling",
// ECCV 2012
// 
// Copyright (c) 2012 Michael Van den Bergh (ETH Zurich). All rights reserved.
// ******************************************************************************

#include "seeds2.h"
#include "math.h"
#include <cstdio>
#include <algorithm>
#include <vector>
#include <iostream>
#include <fstream>
#include "../my_opencv.h"

// Choose color space, LAB is better, HSV is faster (initialization)
//#define LAB_COLORSPACE
//#define HSV_COLORSPACE
#define YUV_COLORSPACE

// Enables means post-processing
#define MEANS

// Enables 3x3 smoothing prior
#define PRIOR

// If enabled, each block-updating level will be iterated twice.
// Disable this to speed up the algorithm.
#define DOUBLE_STEPS
#define REQ_CONF 0.1

#define MINIMUM_NR_SUBLABELS 1



void SEEDS::run() 
{
	// block updates
	while (seeds_current_level >= 0)
	{
		#ifdef DOUBLE_STEPS
			update_blocks(seeds_current_level, REQ_CONF);
		#endif
		update_blocks(seeds_current_level);
		seeds_current_level = go_down_one_level();
	}

	#ifdef MEANS
		compute_means();
		update_pixels_means();
		update_pixels_means();
		update_pixels_means();
		update_pixels_means();
	#else
		update_pixels();
		update_pixels();
		update_pixels();
		update_pixels();
	#endif


	for (int y=0; y<img_height; y++)
	{
		for (int x=0; x<img_width; x++) 
		{
			this->p_s_pix_map[y*img_width+x] = labels[seeds_top_level][y*img_width+x];
		}
	}
	valid_s_pix_num = count_superpixels();

	//statistic_analysis_for_spix(p_spix, valid_s_pix_num);

	//DrawContoursAroundSegments(this->p_edge_map, this->p_s_pix_map, img_width, img_height, 255);  //debug
	////for (int y=0; y<img_height; y++)
	////{
	////	for (int x=0; x<img_width; x++) 
	////	{
	////		if (this->p_edge_map[y * img_width + x] > 0)
	////		{
	////			int index24bit = y * img_width* 3 + x * 3;
	////			this->p_src_img[index24bit + 0] = 255;
	////			this->p_src_img[index24bit + 1] = 128;
	////			this->p_src_img[index24bit + 2] = 128;
	////		}
	////	}
	////}

}


void SEEDS::release_mem()
{
	

	if (initialized)
	{
		for (int level=0; level<seeds_nr_levels; level++)
		{
			//for (int label=0; label<nr_labels[level]; label++)
			//{
			//	delete[] histogram[level][label];
			//}

			delete[] labels[level];
			delete[] parent[level];
			delete[] nr_partitions[level];

		}
		delete[] labels;
		delete[] parent;
		delete[] nr_partitions;
		delete[] nr_labels;
		delete[] nr_w;
		delete[] nr_h;

#ifdef LAB_COLORSPACE
		delete[] bin_cutoff1;
		delete[] bin_cutoff2;
		delete[] bin_cutoff3;
#endif

#ifdef MEANS
		delete[] L_channel;
		delete[] A_channel;
		delete[] B_channel;
#endif

		//delete[] p_spix;
	}
}

void unsigned_char_2_iplimage(unsigned char* p_src_img, IplImage* p_cv_img)
{
	int x = 0, y = 0;
	unsigned char* p_cv_data = (unsigned char*)p_cv_img->imageData;
	int linebyte = p_cv_img->widthStep;
	int index = 0;

	if (p_cv_img->nChannels == 3)
	{
		for (y = 0; y < p_cv_img->height; y++)
		{
			for (x = 0; x < p_cv_img->width; x++)
			{
				index = y * linebyte + x * 3;
				p_cv_data[index + 0] = p_src_img[index + 0];
				p_cv_data[index + 1] = p_src_img[index + 1];
				p_cv_data[index + 2] = p_src_img[index + 2];
			}
		}
	}
	else
	{
		for (y = 0; y < p_cv_img->height; y++)
		{
			for (x = 0; x < p_cv_img->width; x++)
			{
				index = y * linebyte + x;
				p_cv_data[index] = p_src_img[index];
			}
		}
	}
}


//void SEEDS::show_img_by_opencv(unsigned char* p_src_img, int width, int height, int bitcount, char* p_win_name, bool is_yuv)
//{
//	IplImage* p_cv_img = cvCreateImage(cvSize(width, height), 8, bitcount >> 3);
//	unsigned_char_2_iplimage(p_src_img, p_cv_img);
//	if (is_yuv)
//	{
//		cvCvtColor(p_cv_img, p_cv_img, CV_YCrCb2BGR);
//	}
//	cvNamedWindow(p_win_name);
//	cvShowImage(p_win_name, p_cv_img);
//	cvReleaseImage(&p_cv_img);
//}


SEEDS::SEEDS()
{
	initialized = false;

	iteration = 0;
	this->img_width = IMG_MAX_SIZE;
	this->img_height = IMG_MAX_SIZE;
	this->nr_channels = 3;
	this->nr_bins = 5;
	this->merge_threshold = merge_threshold;
	this->min_size = 0;
	this->seeds_nr_levels = 6;
	forwardbackward = true;
	histogram_size = nr_bins*nr_bins*nr_bins;



	image_bins = new UINT[IMG_MAX_SIZE*IMG_MAX_SIZE];
	image_l = new float[IMG_MAX_SIZE*IMG_MAX_SIZE];
	image_a = new float[IMG_MAX_SIZE*IMG_MAX_SIZE];
	image_b = new float[IMG_MAX_SIZE*IMG_MAX_SIZE];
	edge_w = new float[IMG_MAX_SIZE*IMG_MAX_SIZE];
	edge_h = new float[IMG_MAX_SIZE*IMG_MAX_SIZE];
	forwardbackward = true;
	histogram_size = nr_bins*nr_bins*nr_bins;


	//Evaluation
	p_edge_map  = new unsigned char[IMG_MAX_SIZE * IMG_MAX_SIZE];
	p_seg_map   = new unsigned char[IMG_MAX_SIZE * 3 * IMG_MAX_SIZE];
    p_s_pix_map = new short[IMG_MAX_SIZE * IMG_MAX_SIZE];
    p_src_img   = new unsigned char[IMG_MAX_SIZE * 3 * IMG_MAX_SIZE];
	p_spix      = new S_Pix[MAX_S_PIX_NUM];

	histogram = new int**[seeds_nr_levels]; // histograms are kept for each level and each label [level][label][bin]
	T = new int*[seeds_nr_levels]; // block sizes are kept at each level [level][label]
	for (int level=0; level<seeds_nr_levels; level++)
	{
		histogram[level] = new int*[img_width*img_height]; // histograms are kept for each level and each label [level][label][bin]
		T[level] = new int[img_width*img_height]; // block sizes are kept at each level [level][label]
		for (int label=0; label<img_width*img_height; label++)
		{
			histogram[level][label] = new int[histogram_size]; // histogram bins
		}
	}
	this->initialized = true;

}


SEEDS::~SEEDS()
{
	// from constructor
	delete[] image_bins;
	delete[] image_l;
	delete[] image_a;
	delete[] image_b;
	delete[] edge_w;
	delete[] edge_h;
	delete[] p_edge_map;
	delete[] p_seg_map;
	delete[] p_s_pix_map;

	if (initialized)
	{
		for (int level=0; level<seeds_nr_levels; level++)
		{
			for (int label=0; label<nr_labels[level]; label++)
			{
				delete[] histogram[level][label];
			}

			delete[]  histogram[level];
			delete[] T[level];
			delete[] labels[level];
			delete[] parent[level];
			delete[] nr_partitions[level];

		}
		delete[] histogram;
		delete[] T;
		delete[] labels;
		delete[] parent;
		delete[] nr_partitions;
		delete[] nr_labels;
		delete[] nr_w;
		delete[] nr_h;

#ifdef LAB_COLORSPACE
		delete[] bin_cutoff1;
		delete[] bin_cutoff2;
		delete[] bin_cutoff3;
#endif


		#ifdef MEANS
		delete[] L_channel;
		delete[] A_channel;
		delete[] B_channel;
		#endif
	}
}


void SEEDS::init(unsigned char* image, int width, int height, int NR_SUPERPIXELS)
{

	// NOTE: the following values are defined for images from the BSD300 or BSD500 data set.
	// If the input image size differs from 480x320, the following values might no longer be 
	// accurate.
	// For more info on how to select the superpixel sizes, please refer to README.TXT.
	int seed_width = 3; int seed_height = 4; int nr_levels = 4;
	if (width >= height)
	{

		if (NR_SUPERPIXELS >= 600) {seed_width = 2; seed_height = 2; nr_levels = 4;}
		if (NR_SUPERPIXELS == 400) {seed_width = 3; seed_height = 2; nr_levels = 4;}
		if (NR_SUPERPIXELS == 266) {seed_width = 3; seed_height = 3; nr_levels = 4;}
		if (NR_SUPERPIXELS == 200) {seed_width = 3; seed_height = 4; nr_levels = 4;}
		if (NR_SUPERPIXELS == 150) {seed_width = 2; seed_height = 2; nr_levels = 5;}
		if (NR_SUPERPIXELS == 100) {seed_width = 3; seed_height = 2; nr_levels = 5;}
		if (NR_SUPERPIXELS == 50)  {seed_width = 3; seed_height = 4; nr_levels = 5;}
		if (NR_SUPERPIXELS == 25)  {seed_width = 3; seed_height = 2; nr_levels = 6;}
		if (NR_SUPERPIXELS == 17)  {seed_width = 3; seed_height = 3; nr_levels = 6;}
		if (NR_SUPERPIXELS == 12)  {seed_width = 3; seed_height = 4; nr_levels = 6;}
		if (NR_SUPERPIXELS == 9)  {seed_width = 2; seed_height = 2; nr_levels = 7;}
		if (NR_SUPERPIXELS == 6)  {seed_width = 3; seed_height = 2; nr_levels = 7;}
	}
	else
	{
		if (NR_SUPERPIXELS >= 600) {seed_width = 2; seed_height = 2; nr_levels = 4;}
		if (NR_SUPERPIXELS == 400) {seed_width = 2; seed_height = 3; nr_levels = 4;}
		if (NR_SUPERPIXELS == 266) {seed_width = 3; seed_height = 3; nr_levels = 4;}
		if (NR_SUPERPIXELS == 200) {seed_width = 4; seed_height = 3; nr_levels = 4;}
		if (NR_SUPERPIXELS == 150) {seed_width = 2; seed_height = 2; nr_levels = 5;}
		if (NR_SUPERPIXELS == 100) {seed_width = 2; seed_height = 3; nr_levels = 5;}
		if (NR_SUPERPIXELS == 50)  {seed_width = 4; seed_height = 3; nr_levels = 5;}
		if (NR_SUPERPIXELS == 25)  {seed_width = 2; seed_height = 3; nr_levels = 6;}
		if (NR_SUPERPIXELS == 17)  {seed_width = 3; seed_height = 3; nr_levels = 6;}
		if (NR_SUPERPIXELS == 12)  {seed_width = 4; seed_height = 3; nr_levels = 6;}
		if (NR_SUPERPIXELS == 9)  {seed_width = 2; seed_height = 2; nr_levels = 7;}
		if (NR_SUPERPIXELS == 6)  {seed_width = 2; seed_height = 3; nr_levels = 7;}
	}

	iteration = 0;

	this->img_width  = width;
	this->img_height = height;

	this->seeds_w = seed_width;
	this->seeds_h = seed_height;
	this->seeds_nr_levels = nr_levels;
	this->seeds_top_level = nr_levels - 1;
	this->seeds_current_level = nr_levels - 2;


	// assign labels
	assign_labels();
	
	// adaptive histograms
	#ifdef LAB_COLORSPACE
		lab_get_histogram_cutoff_values(image);
	#endif

	// convert input image to LAB or HSV
	for (int x=0; x<width; x++)
		for (int y=0; y<height; y++)
		{
			int i = y*width + x;
			int r = (image[i] >> 16) & 0xFF;
			int g = (image[i] >>  8) & 0xFF;
			int b = (image[i]      ) & 0xFF;
			float L;
			float A;
			float B; 
			#ifdef LAB_COLORSPACE
				image_bins[i] = RGB2LAB_special(r, g, b, &L, &A, &B);
				image_l[i] = L/100.0;
				image_a[i] = (A+128.0)/255.0;
				image_b[i] = (B+128.0)/255.0;
			#endif
			#ifdef HSV_COLORSPACE
				image_bins[i] = RGB2HSV(r, g, b, &L, &A, &B);
                image_l[i] = L;
				image_a[i] = A;
				image_b[i] = B;
			#endif
			#ifdef YUV_COLORSPACE
				int index_24bit = y * width * 3 + x * 3 ;

				float Y = (float)image[index_24bit + 0] / 255.0f;
				float U = (float)image[index_24bit + 1] / 255.0f;
				float V = (float)image[index_24bit + 2] / 255.0f;

				image_bins[i] =  (Y * nr_bins) + nr_bins*(U * nr_bins) + nr_bins*nr_bins*(V * nr_bins);
				image_l[i] = Y;
				image_a[i] = U;
				image_b[i] = V;

				p_src_img[index_24bit + 0]   = image[index_24bit + 0];
				p_src_img[index_24bit + 1]   = image[index_24bit + 1];
				p_src_img[index_24bit + 2]   = image[index_24bit + 2];

			#endif
		}

	// compute histograms
	compute_histograms();



}




// Assign Labels
// adds labeling to all the blocks at all levels and sets the correct parents
void SEEDS::assign_labels()
{
	// initialize lists
	labels = new UINT*[seeds_nr_levels];
	parent = new UINT*[seeds_nr_levels];
	nr_partitions = new UINT*[seeds_nr_levels];
	nr_labels = new UINT[seeds_nr_levels];
	nr_w = new int[seeds_nr_levels];
	nr_h = new int[seeds_nr_levels];

	// level 0 (base level): seeds_w x seeds_h
	int level = 0;
	int nr_seeds_w = floor(img_width/seeds_w);
	int nr_seeds_h = floor(img_height/seeds_h);
	int step_w = seeds_w;
	int step_h = seeds_h;
	int nr_seeds = nr_seeds_w*nr_seeds_h;
	nr_labels[level] = nr_seeds;
	nr_w[level] = nr_seeds_w;
	nr_h[level] = nr_seeds_h;
	labels[level] = new UINT[img_width*img_height];
	parent[level] = new UINT[nr_seeds];
	nr_partitions[level] = new UINT[nr_seeds];
	for (int i=0; i<nr_seeds; i++) nr_partitions[level][i] = 1; // base level: no further partitioning
	for (int y=0; y<img_height; y++)
		for (int x=0; x<img_width; x++)
		{
			int label_x = floor((float)x/step_w);
			if (label_x >= nr_seeds_w) label_x = nr_seeds_w-1;
			int label_y = floor((float)y/step_h);
			if (label_y >= nr_seeds_h) label_y = nr_seeds_h-1;
			labels[level][y*img_width+x] = label_y*nr_seeds_w + label_x;
			// parents will be set in the next loop
		}

	// level 1, 2, 3 ... until top level
	for (int level = 1; level < seeds_nr_levels; level++)
	{
		nr_seeds_w /= 2; // always partitioned in 2x2 sub-blocks
		nr_seeds_h /= 2; // always partitioned in 2x2 sub-blocks
		step_w *= 2;
		step_h *= 2;
		nr_seeds = nr_seeds_w*nr_seeds_h;
		nr_labels[level] = nr_seeds;
		nr_w[level] = nr_seeds_w;
		nr_h[level] = nr_seeds_h;
		labels[level] = new UINT[img_width*img_height];
		parent[level] = new UINT[nr_seeds];
		nr_partitions[level] = new UINT[nr_seeds];
		for (int i=0; i<nr_seeds; i++) nr_partitions[level][i] = 0; // will be added in add_blocks
		for (int y=0; y<img_height; y++)
			for (int x=0; x<img_width; x++)
			{
				int label_x = floor((float)x/step_w);
				if (label_x >= nr_seeds_w) label_x = nr_seeds_w-1;
				int label_y = floor((float)y/step_h);
				if (label_y >= nr_seeds_h) label_y = nr_seeds_h-1;
				labels[level][y*img_width+x] = label_y*nr_seeds_w + label_x;
				parent[level-1][labels[level-1][y*img_width+x]] = labels[level][y*img_width+x]; // set parent
			}
	}

	seeds_top_level = seeds_nr_levels - 1;
}


//=================================================================================
/// DrawContoursAroundSegments
///
/// Internal contour drawing option exists. One only needs to comment the if
/// statement inside the loop that looks at neighbourhood.
//=================================================================================
void SEEDS::DrawContoursAroundSegments(
	unsigned char*			ubuff,
	const short*			labels,
	const int&				width,
	const int&				height,
	const unsigned int&		color )
{
	const int dx8[8] = {-1, -1,  0,  1, 1, 1, 0, -1};
	const int dy8[8] = { 0, -1, -1, -1, 0, 1, 1,  1};

	int sz = width*height;

	vector<bool> istaken(sz, false);

	int mainindex(0);
	for( int j = 0; j < height; j++ )
	{
		for( int k = 0; k < width; k++ )
		{

			int np(0);
			ubuff[mainindex] = 0;
			for( int i = 0; i < 8; i++ )
			{
				int x = k + dx8[i];
				int y = j + dy8[i];

				if( (x >= 0 && x < width) && (y >= 0 && y < height) )
				{
					int index = y*width + x;

					if( false == istaken[index] )//comment this to obtain internal contours
					{
						if( labels[mainindex] != labels[index] ) np++;
					}
				}
			}

			if( np > 1 )//change to 2 or 3 for thinner lines
			{
				ubuff[mainindex] = 255;
				istaken[mainindex] = true;
			}
			mainindex++;

		}
	}
}

void SEEDS::lab_get_histogram_cutoff_values(UINT* image)
{
	// get image lists and histogram cutoff values
	vector<float> list_channel1;
	vector<float> list_channel2;
	vector<float> list_channel3;
	vector<float>::iterator it;
	int samp = 5;
	int ctr = 0;
	for (int x=0; x<img_width; x+=samp)
		for (int y=0; y<img_height; y+=samp)
		{
			int i = y*img_width +x;
			int r = (image[i] >> 16) & 0xFF;
			int g = (image[i] >>  8) & 0xFF;
			int b = (image[i]      ) & 0xFF;
			float L = 0;
			float A = 0;
			float B = 0;
			RGB2LAB(r, g, b, &L, &A, &B);
			list_channel1.push_back(L);
			list_channel2.push_back(A);
			list_channel3.push_back(B);
			ctr++;
		}
	bin_cutoff1 = new float[nr_bins];
	bin_cutoff2 = new float[nr_bins];
	bin_cutoff3 = new float[nr_bins];
	for (int i=1; i<nr_bins; i++)
	{
		int N = (int) floor((float) (i*ctr)/ (float)nr_bins);
		nth_element(list_channel1.begin(), list_channel1.begin()+N, list_channel1.end());
		it = list_channel1.begin()+N;
		bin_cutoff1[i-1] = *it;
		nth_element(list_channel2.begin(), list_channel2.begin()+N, list_channel2.end());
		it = list_channel2.begin()+N;
		bin_cutoff2[i-1] = *it;
		nth_element(list_channel3.begin(), list_channel3.begin()+N, list_channel3.end());
		it = list_channel3.begin()+N;
		bin_cutoff3[i-1] = *it;
	}
	bin_cutoff1[nr_bins-1] = 300;
	bin_cutoff2[nr_bins-1] = 300;
	bin_cutoff3[nr_bins-1] = 300;
}


void SEEDS::compute_means()
{
	L_channel = new float[nr_labels[seeds_top_level]];
	A_channel = new float[nr_labels[seeds_top_level]];
	B_channel = new float[nr_labels[seeds_top_level]];

	// clear counted LAB values
	for (int label=0; label<nr_labels[seeds_top_level]; label++)
	{
		L_channel[label] = 0;
		A_channel[label] = 0;
		B_channel[label] = 0;
	}

	// sweep image
	for (int i=0; i<img_width*img_height; i++)
	{
		int label = labels[seeds_top_level][i];
		// add LAB values to total
		L_channel[label] += image_l[i];
		A_channel[label] += image_a[i];
		B_channel[label] += image_b[i];
	}
}

void SEEDS::compute_histograms(int until_level)
{
	if (until_level == -1) until_level = seeds_nr_levels - 1;
	until_level++;

	// initialize histograms in the first iteration
	if (iteration == 0)
	{
		//histogram = new int**[seeds_nr_levels]; // histograms are kept for each level and each label [level][label][bin]
		//T = new int*[seeds_nr_levels]; // block sizes are kept at each level [level][label]
		//for (int level=0; level<seeds_nr_levels; level++)
		//{
		//	histogram[level] = new int*[nr_labels[level]]; // histograms are kept for each level and each label [level][label][bin]
		//	T[level] = new int[nr_labels[level]]; // block sizes are kept at each level [level][label]
		//	for (int label=0; label<nr_labels[level]; label++)
		//	{
		//		histogram[level][label] = new int[histogram_size]; // histogram bins
		//	}
		//}
	}

	// clear histograms
	for (int level=0; level<seeds_nr_levels; level++)
		for (int label=0; label<nr_labels[level]; label++)
		{
			for (int j=0; j<histogram_size; j++)
				histogram[level][label][j] = 0;
			T[level][label] = 0;
		}

	// build histograms on the first level by adding the pixels to the blocks
	for (int x=0; x<img_width; x++)
		for (int y=0; y<img_height; y++)
		{					
			int i = y*img_width +x;
			add_pixel(0, labels[0][i], x, y);
		}

	// build histograms on the upper levels by adding the histogram from the level below
	for (int level=1; level<until_level; level++)
	{
		for (int label=0; label<nr_labels[level-1]; label++)
		{
			add_block(level, parent[level-1][label], level-1, label);
		}
	}
}

void SEEDS::compute_histograms_ex()
{
	// clear histograms
	for (int level=0; level<seeds_nr_levels; level++)
		for (int label=0; label<nr_labels[level]; label++)
		{
			for (int j=0; j<histogram_size; j++)
				histogram[level][label][j] = 0;
			T[level][label] = 0;
		}

	for (int level=0; level<seeds_nr_levels; level++)
		for (int x=0; x<img_width; x++)
			for (int y=0; y<img_height; y++)
			{					
				int i = y*img_width +x;
				//add_pixel(level, labels[level][i], x, y);
				histogram[level][labels[level][i]][image_bins[y*img_width+x]]++;
				T[level][labels[level][i]]++;
			}

}

void SEEDS::compute_edges()
{
	// compute edges
	/*for (int x=0;x<img_width-1; x++)
		for (int y=0; y<height; y++)
		{
			float Ex = image_l[y*img_width+x] - image_l[y*img_width+(x+1)];
			#ifdef HSV_COLORSPACE
				if (Ex > 0.5) Ex = 1.0 - Ex;
			#endif
			float Elx = image_a[y*img_width+x] - image_a[y*img_width+(x+1)];
			float Ellx = image_b[y*img_width+x] - image_b[y*img_width+(x+1)];
			edge_w[y*img_width+x] = sqrt(Ex*Ex + Elx*Elx + Ellx*Ellx);
		}
	for (int y=0;y<height-1; y++)
		for (int x=0; x<img_width; x++)
		{
			float Ey = image_l[y*img_width+x] - image_l[(y+1)*img_width+x];
			#ifdef HSV_COLORSPACE
				if (Ey > 0.5) Ey = 1.0 - Ey;
			#endif
			float Ely = image_a[y*img_width+x] - image_a[(y+1)*img_width+x];
			float Elly = image_b[y*img_width+x] - image_b[(y+1)*img_width+x];
			edge_h[y*img_width+x] = sqrt(Ey*Ey + Ely*Ely + Elly*Elly);
		}*/
}



// Merging Algorithm
void SEEDS::update_blocks(int level, float req_confidence)
{
	int labelA;
	int labelB;
	int sublabel;
	int priorA;
	int priorB;
	bool done;

	int step = nr_w[level];
	
	// horizontal bidirectional block updating
	for (int y=1; y<nr_h[level]-2; y++)
		for (int x=1; x<nr_w[level]-2; x++) 
		{
			// choose a label at the current level
			sublabel = y*step+x; 
			// get the label at the top level (= superpixel label)
			labelA = parent[level][y*step+x]; 
			// get the neighboring label at the top level (= superpixel label)
			labelB = parent[level][y*step+x+1]; 

			if (labelA != labelB)
			{
				// get the surrounding labels at the top level, to check for splitting 
				int a11 = parent[level][(y-1)*step+(x-1)]; 
				int a12 = parent[level][(y-1)*step+(x)];
				int a13 = parent[level][(y-1)*step+(x+1)];
				int a14 = parent[level][(y-1)*step+(x+2)];
				int a21 = parent[level][(y)*step+(x-1)];
				int a22 = parent[level][(y)*step+(x)];
				int a23 = parent[level][(y)*step+(x+1)];
				int a24 = parent[level][(y)*step+(x+2)];
				int a31 = parent[level][(y+1)*step+(x-1)];
				int a32 = parent[level][(y+1)*step+(x)];
				int a33 = parent[level][(y+1)*step+(x+1)];
				int a34 = parent[level][(y+1)*step+(x+2)];

				done = false;

				if (nr_partitions[seeds_top_level][labelA] > MINIMUM_NR_SUBLABELS)
				{
					// remove sublabel from A, get intersections with A and B
					if (nr_partitions[seeds_top_level][labelA] <= 2) // == 2
					{
						// run algorithm as usual
						delete_block(seeds_top_level, labelA, level, sublabel);
						float intA = intersection(seeds_top_level, labelA, level, sublabel);
						float intB = intersection(seeds_top_level, labelB, level, sublabel);
						float confidence = fabs(intA - intB);
						// add to label with highest intersection
						if ((intB > intA) && (confidence > req_confidence))
						{
							add_block(seeds_top_level, labelB, level, sublabel);
							done = true;
						}
						else
						{
							add_block(seeds_top_level, labelA, level, sublabel);
						}
					}
					else if (nr_partitions[seeds_top_level][labelA] > 2) // 3 or more partitions
					{
						if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, true, true))
						{
							// run algorithm as usual
							delete_block(seeds_top_level, labelA, level, sublabel);
							float intA = intersection(seeds_top_level, labelA, level, sublabel);
							float intB = intersection(seeds_top_level, labelB, level, sublabel);
							float confidence = fabs(intA - intB);
							// add to label with highest intersection
							if ((intB > intA) && (confidence > req_confidence))
							{
								add_block(seeds_top_level, labelB, level, sublabel);
								done = true;
							}
							else
							{
								add_block(seeds_top_level, labelA, level, sublabel);
							}
						}
					}
				}

				if ((!done) && (nr_partitions[seeds_top_level][labelB] > MINIMUM_NR_SUBLABELS))
				{
					// try opposite direction
					sublabel = y*step+x+1;
					if (nr_partitions[seeds_top_level][labelB] <= 2) // == 2
					{
						// run algorithm as usual
						delete_block(seeds_top_level, labelB, level, sublabel);
						float intA = intersection(seeds_top_level, labelA, level, sublabel);
						float intB = intersection(seeds_top_level, labelB, level, sublabel);
						float confidence = fabs(intA - intB);
						if ((intA > intB) && (confidence > req_confidence))
						{
							add_block(seeds_top_level, labelA, level, sublabel);
							x++;
						}
						else
						{
							add_block(seeds_top_level, labelB, level, sublabel);
						}
					}
					else if (nr_partitions[seeds_top_level][labelB] > 2)
					{
						if (!check_split(a12, a13, a14, a22, a23, a24, a32, a33, a34, true, false))
						{
							// run algorithm as usual
							delete_block(seeds_top_level, labelB, level, sublabel);
							float intA = intersection(seeds_top_level, labelA, level, sublabel);
							float intB = intersection(seeds_top_level, labelB, level, sublabel);
							float confidence = fabs(intA - intB);
							if ((intA > intB) && (confidence > req_confidence))
							{
								add_block(seeds_top_level, labelA, level, sublabel);
								x++;
							}
							else
							{
								add_block(seeds_top_level, labelB, level, sublabel);
							}
						}
					}
				}											
			}
		}
		update_labels(level);

	// vertical bidirectional
	for (int x=1; x<nr_w[level]-2; x++)
		for (int y=1; y<nr_h[level]-2; y++)
		{
			// choose a label at the current level
			sublabel = y*step+x; 
			// get the label at the top level (= superpixel label)
			labelA = parent[level][y*step+x]; 
			// get the neighboring label at the top level (= superpixel label)
			labelB = parent[level][(y+1)*step+x]; 

			if (labelA != labelB)
			{
				int a11 = parent[level][(y-1)*step+(x-1)];
				int a12 = parent[level][(y-1)*step+(x)];
				int a13 = parent[level][(y-1)*step+(x+1)];
				int a21 = parent[level][(y)*step+(x-1)];
				int a22 = parent[level][(y)*step+(x)];
				int a23 = parent[level][(y)*step+(x+1)];
				int a31 = parent[level][(y+1)*step+(x-1)];
				int a32 = parent[level][(y+1)*step+(x)];
				int a33 = parent[level][(y+1)*step+(x+1)];
				int a41 = parent[level][(y+2)*step+(x-1)];
				int a42 = parent[level][(y+2)*step+(x)];
				int a43 = parent[level][(y+2)*step+(x+1)];

				done = false;
				if (nr_partitions[seeds_top_level][labelA] > MINIMUM_NR_SUBLABELS)
				{
					// remove sublabel from A, get intersections with A and B
					if (nr_partitions[seeds_top_level][labelA] <= 2) // == 2
					{
						// run algorithm as usual
						delete_block(seeds_top_level, labelA, level, sublabel);
						float intA = intersection(seeds_top_level, labelA, level, sublabel);
						float intB = intersection(seeds_top_level, labelB, level, sublabel);
						float confidence = fabs(intA - intB);
						// add to label with highest intersection
						if ((intB > intA) && (confidence > req_confidence))
						{
							add_block(seeds_top_level, labelB, level, sublabel);
							//y++;
							done = true;
						}
						else
						{
							add_block(seeds_top_level, labelA, level, sublabel);
						}
					}
					else if (nr_partitions[seeds_top_level][labelA] > 2) // 3 or more partitions
					{
						if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, false, true))
						{
							// run algorithm as usual
							delete_block(seeds_top_level, labelA, level, sublabel);
							float intA = intersection(seeds_top_level, labelA, level, sublabel);
							float intB = intersection(seeds_top_level, labelB, level, sublabel);
							float confidence = fabs(intA - intB);
							// add to label with highest intersection
							if ((intB > intA) && (confidence > req_confidence))
							{
								add_block(seeds_top_level, labelB, level, sublabel);
								//y++;
								done = true;
							}
							else
							{
								add_block(seeds_top_level, labelA, level, sublabel);
							}
						}
					}
				}

				if ((!done) && (nr_partitions[seeds_top_level][labelB] > MINIMUM_NR_SUBLABELS))
				{
					// try opposite direction
					sublabel = (y+1)*step+x;
					if (nr_partitions[seeds_top_level][labelB] <= 2) // == 2
					{
						// run algorithm as usual
						delete_block(seeds_top_level, labelB, level, sublabel);
						float intA = intersection(seeds_top_level, labelA, level, sublabel);
						float intB = intersection(seeds_top_level, labelB, level, sublabel);
						float confidence = fabs(intA - intB);
						if ((intA > intB) && (confidence > req_confidence))
						{
							add_block(seeds_top_level, labelA, level, sublabel);
							y++;
						}
						else
						{
							add_block(seeds_top_level, labelB, level, sublabel);
						}
					}
					else if (nr_partitions[seeds_top_level][labelB] > 2)
					{
						if (!check_split(a21, a22, a23, a31, a32, a33, a41, a42, a43, false, false))
						{
							// run algorithm as usual
							delete_block(seeds_top_level, labelB, level, sublabel);
							float intA = intersection(seeds_top_level, labelA, level, sublabel);
							float intB = intersection(seeds_top_level, labelB, level, sublabel);
							float confidence = fabs(intA - intB);
							if ((intA > intB) && (confidence > req_confidence))
							{
								add_block(seeds_top_level, labelA, level, sublabel);
								y++;
							}
							else
							{
								add_block(seeds_top_level, labelB, level, sublabel);
							}
						}
					}
				}											
			}
		}
		update_labels(level);
}



int SEEDS::go_down_one_level()
{
	int old_level = seeds_current_level;
	int new_level = seeds_current_level - 1;

	if (new_level < 0) return -1;

	// go through labels of top level
	for (int x=0; x<nr_w[seeds_top_level]; x++)
		for (int y=0; y<nr_h[seeds_top_level]; y++)
		{
			// reset nr_partitions
			nr_partitions[seeds_top_level][y*nr_w[seeds_top_level]+x] = 0;
		}

	// go through labels of new_level
	for (int x=0; x<nr_w[new_level]; x++)
		for (int y=0; y<nr_h[new_level]; y++)
		{
			// assign parent = parent of old_label
			int p = parent[old_level][parent[new_level][y*nr_w[new_level]+x]];
			parent[new_level][y*nr_w[new_level]+x] = p;
			// add nr_partitions[label] to parent
			nr_partitions[seeds_top_level][p] += nr_partitions[new_level][y*nr_w[new_level]+x];
		}

	return new_level;
}



// Border Updating Algorithm
void SEEDS::update_pixels()
{
	int labelA, labelA1, labelA2, labelA3;
	int labelB, labelB1, labelB2, labelB3;
	int priorA = 0;
	int priorB = 0;
	float edgeA = 0;
	float edgeB = 0;

	if (forwardbackward)
	{
		forwardbackward = false;
		for (int y=1; y<img_height-1; y++)
			for (int x=1; x<img_width-1; x++) 
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a14 = labels[seeds_top_level][(y-1)*img_width+(x+2)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a24 = labels[seeds_top_level][(y)*img_width+(x+2)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)];
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a34 = labels[seeds_top_level][(y+1)*img_width+(x+2)];

				// horizontal bidirectional
				labelA = a22;
				labelB = a23;
				if (labelA != labelB)
				{
					if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, true, true))
					{
						#ifdef PRIOR
						priorA = threebyfour(x,y,labelA);
						priorB = threebyfour(x,y,labelB);
						#endif

						if (probability(image_bins[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
						{
							update(seeds_top_level, labelB, x, y);
						}
						else if (!check_split(a12, a13, a14, a22, a23, a24, a32, a33, a34, true, false))
						{
							if (probability(image_bins[y*img_width+x+1], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
							{
								update(seeds_top_level, labelA, x+1, y);
								x++;
							}
						}
					}
				}
			}


		for (int x=1; x<img_width-1; x++)
			for (int y=1; y<img_height-2; y++)
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)]; 
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a41 = labels[seeds_top_level][(y+2)*img_width+(x-1)];
				int a42 = labels[seeds_top_level][(y+2)*img_width+(x)];
				int a43 = labels[seeds_top_level][(y+2)*img_width+(x+1)];

				// vertical bidirectional
				labelA = a22;
				labelB = a32;
				if (labelA != labelB)
				{
					if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, false, true))
					{
						#ifdef PRIOR
						priorA = fourbythree(x,y,labelA);
						priorB =  fourbythree(x,y,labelB);
						#endif

						if (probability(image_bins[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
						{
							update(seeds_top_level, labelB, x, y);
						} 
						else if (!check_split(a21, a22, a23, a31, a32, a33, a41, a42, a43, false, false))
						{
							if (probability(image_bins[(y+1)*img_width+x], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
							{
								update(seeds_top_level, labelA, x, y+1);
								y++;
							}
						}
					}
				}
			}
	}
	else
	{
		forwardbackward = true;
		for (int y=1; y<img_height-1; y++)
			for (int x=1; x<img_width-2; x++) 
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a14 = labels[seeds_top_level][(y-1)*img_width+(x+2)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a24 = labels[seeds_top_level][(y)*img_width+(x+2)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)];
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a34 = labels[seeds_top_level][(y+1)*img_width+(x+2)];

				// horizontal bidirectional
				labelA = a22;
				labelB = a23;
				if (labelA != labelB)
				{
					if (!check_split(a12, a13, a14, a22, a23, a24, a32, a33, a34, true, false))
					{
						#ifdef PRIOR
						priorA = threebyfour(x,y,labelA);
						priorB = threebyfour(x,y,labelB);
						#endif

						if (probability(image_bins[y*img_width+x+1], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
						{
							update(seeds_top_level, labelA, x+1, y);
							x++;
						}
						else if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, true, true))
						{
							if (probability(image_bins[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
							{
								update(seeds_top_level, labelB, x, y);
							}
						}
					}
				}
			}


		for (int x=1; x<img_width-1; x++)
			for (int y=1; y<img_height-2; y++)
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)]; 
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a41 = labels[seeds_top_level][(y+2)*img_width+(x-1)];
				int a42 = labels[seeds_top_level][(y+2)*img_width+(x)];
				int a43 = labels[seeds_top_level][(y+2)*img_width+(x+1)];

				// vertical bidirectional
				labelA = a22;
				labelB = a32;
				if (labelA != labelB)
				{
					if (!check_split(a21, a22, a23, a31, a32, a33, a41, a42, a43, false, false))
					{
						#ifdef PRIOR						
						priorA = fourbythree(x,y,labelA);
						priorB =  fourbythree(x,y,labelB);
						#endif

						if (probability(image_bins[(y+1)*img_width+x], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
						{
							update(seeds_top_level, labelA, x, y+1);
							y++;
						}
						else if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, false, true))
						{
							if (probability(image_bins[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
							{
								update(seeds_top_level, labelB, x, y);
							}
						}
					}
				}
			}

	}

	// update border pixels
	for (int x=0; x<img_width; x++)
	{
		labelA = labels[seeds_top_level][x];
		labelB = labels[seeds_top_level][img_width+x];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, x, 0);	
		}
		labelA = labels[seeds_top_level][(img_height-1)*img_width + x];
		labelB = labels[seeds_top_level][(img_height-2)*img_width + x];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, x, img_height-1);			
		}
	}
	for (int y=0; y<img_height; y++)
	{
		labelA = labels[seeds_top_level][y*img_width];
		labelB = labels[seeds_top_level][y*img_width+1];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, 0, y);	
		}
		labelA = labels[seeds_top_level][y*img_width + img_width - 1];
		labelB = labels[seeds_top_level][y*img_width + img_width - 2];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, img_width-1, y);			
		}
	}
}


// Border Updating Algorithm
void SEEDS::update_pixels_means()
{
	int labelA = 0, labelA1 = 0, labelA2 = 0, labelA3 = 0;
	int labelB = 0, labelB1 = 0, labelB2 = 0, labelB3 = 0;
	int priorA = 0;
	int priorB = 0;
	float edgeA = 0;
	float edgeB = 0;

	if (forwardbackward)
	{
		forwardbackward = false;
		for (int y=1; y<img_height-1; y++)
			for (int x=1; x<img_width-2; x++) 
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a14 = labels[seeds_top_level][(y-1)*img_width+(x+2)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)];
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)]; 
				int a24 = labels[seeds_top_level][(y)*img_width+(x+2)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)];
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a34 = labels[seeds_top_level][(y+1)*img_width+(x+2)];

				// horizontal bidirectional
				labelA = a22;
				labelB = a23;
				if (labelA != labelB)
				{
					if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, true, true))
					{
						#ifdef PRIOR
						priorA = threebyfour(x,y,labelA);
						priorB = threebyfour(x,y,labelB);
						#endif

						if (probability_means(image_l[y*img_width+x], image_a[y*img_width+x], image_b[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
						{
							update(seeds_top_level, labelB, x, y);
						}
						else if (!check_split(a12, a13, a14, a22, a23, a24, a32, a33, a34, true, false))
						{
							if (probability_means(image_l[y*img_width+x+1], image_a[y*img_width+x+1], image_b[y*img_width+x+1], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
							{
								update(seeds_top_level, labelA, x+1, y);
								x++;
							}
						}
					}
				}
			}


		for (int x=1; x<img_width-1; x++)
			for (int y=1; y<img_height-2; y++)
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)]; 
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a41 = labels[seeds_top_level][(y+2)*img_width+(x-1)];
				int a42 = labels[seeds_top_level][(y+2)*img_width+(x)];
				int a43 = labels[seeds_top_level][(y+2)*img_width+(x+1)];

				// vertical bidirectional
				labelA = a22;
				labelB = a32;
				if (labelA != labelB)
				{
					if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, false, true))
					{
						#ifdef PRIOR
						priorA = fourbythree(x,y,labelA);
						priorB =  fourbythree(x,y,labelB);
						#endif

						if (probability_means(image_l[y*img_width+x], image_a[y*img_width+x], image_b[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
						{
							update(seeds_top_level, labelB, x, y);
						} 
						else if (!check_split(a21, a22, a23, a31, a32, a33, a41, a42, a43, false, false))
						{
							if (probability_means(image_l[(y+1)*img_width+x], image_a[(y+1)*img_width+x], image_b[(y+1)*img_width+x], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
							{
								update(seeds_top_level, labelA, x, y+1);
								y++;
							}
						}
					}
				}
			}
	}
	else
	{
		forwardbackward = true;
		for (int y=1; y<img_height-1; y++)
			for (int x=1; x<img_width-2; x++) 
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a14 = labels[seeds_top_level][(y-1)*img_width+(x+2)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)]; 
				int a24 = labels[seeds_top_level][(y)*img_width+(x+2)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)];
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a34 = labels[seeds_top_level][(y+1)*img_width+(x+2)];

				// horizontal bidirectional
				labelA = a22;
				labelB = a23;
				if (labelA != labelB)
				{
					if (!check_split(a12, a13, a14, a22, a23, a24, a32, a33, a34, true, false))
					{
						#ifdef PRIOR
						priorA = threebyfour(x,y,labelA);
						priorB = threebyfour(x,y,labelB);
						#endif

						if (probability_means(image_l[y*img_width+x+1], image_a[y*img_width+x+1], image_b[y*img_width+x+1], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
						{
							update(seeds_top_level, labelA, x+1, y);
							x++;
						}
						else if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, true, true))
						{
							if (probability_means(image_l[y*img_width+x], image_a[y*img_width+x], image_b[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
							{
								update(seeds_top_level, labelB, x, y);
							}
						}
					}
				}
			}


		for (int x=1; x<img_width-1; x++)
			for (int y=1; y<img_height-2; y++)
			{
				int a11 = labels[seeds_top_level][(y-1)*img_width+(x-1)];
				int a12 = labels[seeds_top_level][(y-1)*img_width+(x)];
				int a13 = labels[seeds_top_level][(y-1)*img_width+(x+1)];
				int a21 = labels[seeds_top_level][(y)*img_width+(x-1)];
				int a22 = labels[seeds_top_level][(y)*img_width+(x)]; 
				int a23 = labels[seeds_top_level][(y)*img_width+(x+1)];
				int a31 = labels[seeds_top_level][(y+1)*img_width+(x-1)];
				int a32 = labels[seeds_top_level][(y+1)*img_width+(x)]; 
				int a33 = labels[seeds_top_level][(y+1)*img_width+(x+1)];
				int a41 = labels[seeds_top_level][(y+2)*img_width+(x-1)];
				int a42 = labels[seeds_top_level][(y+2)*img_width+(x)];
				int a43 = labels[seeds_top_level][(y+2)*img_width+(x+1)];

				// vertical bidirectional
				labelA = a22;
				labelB = a32;
				if (labelA != labelB)
				{
					if (!check_split(a21, a22, a23, a31, a32, a33, a41, a42, a43, false, false))
					{
						#ifdef PRIOR
						priorA = fourbythree(x,y,labelA);
						priorB =  fourbythree(x,y,labelB);
						#endif

						if (probability_means(image_l[(y+1)*img_width+x], image_a[(y+1)*img_width+x], image_b[(y+1)*img_width+x], labelB, labelA, priorB, priorA, edgeB, edgeA)) 
						{
							update(seeds_top_level, labelA, x, y+1);
							y++;
						}
						else if (!check_split(a11, a12, a13, a21, a22, a23, a31, a32, a33, false, true))
						{
							if (probability_means(image_l[y*img_width+x], image_a[y*img_width+x], image_b[y*img_width+x], labelA, labelB, priorA, priorB, edgeA, edgeB)) 
							{
								update(seeds_top_level, labelB, x, y);
							}
						}
					}
				}
			}

	}

	// update border pixels
	for (int x=0; x<img_width; x++)
	{
		labelA = labels[seeds_top_level][x];
		labelB = labels[seeds_top_level][img_width+x];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, x, 0);	
		}
		labelA = labels[seeds_top_level][(img_height-1)*img_width + x];
		labelB = labels[seeds_top_level][(img_height-2)*img_width + x];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, x, img_height-1);			
		}
	}
	for (int y=0; y<img_height; y++)
	{
		labelA = labels[seeds_top_level][y*img_width];
		labelB = labels[seeds_top_level][y*img_width+1];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, 0, y);	
		}
		labelA = labels[seeds_top_level][y*img_width + img_width - 1];
		labelB = labels[seeds_top_level][y*img_width + img_width - 2];
		if (labelA != labelB)
		{
				update(seeds_top_level, labelB, img_width-1, y);			
		}
	}
}



void SEEDS::update(int level, int label_new, int x, int y)
{
	int label_old = labels[level][y*img_width+x];
    delete_pixel_m(level, label_old, x, y);
    add_pixel_m(level, label_new, x, y);
	labels[level][y*img_width+x] = label_new;
}

void SEEDS::add_pixel(int level, int label, int x, int y)
{
	histogram[level][label][image_bins[y*img_width+x]]++;
	T[level][label]++;
}

void SEEDS::add_pixel_m(int level, int label, int x, int y)
{
	histogram[level][label][image_bins[y*img_width+x]]++;
	T[level][label]++;

	#ifdef MEANS
		L_channel[label] += image_l[y*img_width + x];
		A_channel[label] += image_a[y*img_width + x];
		B_channel[label] += image_b[y*img_width + x];
    #endif
}

void SEEDS::delete_pixel(int level, int label, int x, int y)
{
	histogram[level][label][image_bins[y*img_width+x]]--;
	T[level][label]--;
}

void SEEDS::delete_pixel_m(int level, int label, int x, int y)
{
	histogram[level][label][image_bins[y*img_width+x]]--;
	T[level][label]--;
	
	#ifdef MEANS
		L_channel[label] -= image_l[y*img_width + x];
		A_channel[label] -= image_a[y*img_width + x];
		B_channel[label] -= image_b[y*img_width + x];
    #endif
}


void SEEDS::add_block(int level, int label, int sublevel, int sublabel)
{
	parent[sublevel][sublabel] = label;

	for (int n=0; n<histogram_size; n++)
	{
		histogram[level][label][n] += histogram[sublevel][sublabel][n];
	}
	T[level][label] += T[sublevel][sublabel];

	nr_partitions[level][label]++;
}



void SEEDS::delete_block(int level, int label, int sublevel, int sublabel)
{
	parent[sublevel][sublabel] = -1;

	for (int n=0; n<histogram_size; n++)
	{
		histogram[level][label][n] -= histogram[sublevel][sublabel][n];
	}
	T[level][label] -= T[sublevel][sublabel];

	nr_partitions[level][label]--;
}



void SEEDS::update_labels(int level)
{
	for (int i=0; i<img_width*img_height; i++)
	{
		labels[seeds_top_level][i] = parent[level][labels[level][i]];
	}
}




bool SEEDS::probability(int color, int label1, int label2, int prior1, int prior2, float edge1, float edge2)
{
	float P_label1 = (float)histogram[seeds_top_level][label1][color] / (float)T[seeds_top_level][label1];
	float P_label2 = (float)histogram[seeds_top_level][label2][color] / (float)T[seeds_top_level][label2];

	#ifdef PRIOR
		P_label1 *= (float) prior1;
		P_label2 *= (float) prior2;
	#else
		P_label1 = (float)histogram[seeds_top_level][label1][color] / (float)T[seeds_top_level][label1];
		P_label2 = (float)histogram[seeds_top_level][label2][color] / (float)T[seeds_top_level][label2];
	#endif

	return (P_label2 > P_label1);
}

bool SEEDS::probability_means(float L, float a, float b, int label1, int label2, int prior1, int prior2, float edge1, float edge2)
{
	#ifdef MEANS
		float L1 = L_channel[label1] / T[seeds_top_level][label1];
		float a1 = A_channel[label1] / T[seeds_top_level][label1];
		float b1 = B_channel[label1] / T[seeds_top_level][label1];
		float L2 = L_channel[label2] / T[seeds_top_level][label2];
		float a2 = A_channel[label2] / T[seeds_top_level][label2];
		float b2 = B_channel[label2] / T[seeds_top_level][label2];

		float P_label1 = (L-L1)*(L-L1) + (a-a1)*(a-a1) + (b-b1)*(b-b1);
		float P_label2 = (L-L2)*(L-L2) + (a-a2)*(a-a2) + (b-b2)*(b-b2);

		#ifdef PRIOR
			P_label1 /= prior1;
			P_label2 /= prior2;
		#endif

		return (P_label1 > P_label2);
	#endif

	return false;
}



void SEEDS::LAB2RGB(float L, float a, float b, int* R, int* G, int* B)
{
	float T1 = 0.008856;
	float T2 = 0.206893;

	float fY = pow((L + 16.0) / 116.0, 3);
	bool YT = fY > T1;
	fY = (!YT) * (L / 903.3) + YT * fY;
	float Y = fY;

	// Alter fY slightly for further calculations
	fY = YT * pow(fY,1.0/3.0) + (!YT) * (7.787 * fY + 16.0/116.0);

	float fX = a / 500.0 + fY;
	bool XT = fX > T2;
	float X = (XT * pow(fX, 3) + (!XT) * ((fX - 16.0/116.0) / 7.787));

	float fZ = fY - b / 200.0;
	bool ZT = fZ > T2;
	float Z = (ZT * pow(fZ, 3) + (!ZT) * ((fZ - 16.0/116.0) / 7.787));

	X = X * 0.950456 * 255.0;
	Y = Y * 255.0;
	Z = Z * 1.088754 * 255.0;

	*R = (int) (3.240479*X - 1.537150*Y - 0.498535*Z);
	*G = (int) (-0.969256*X + 1.875992*Y + 0.041556*Z);
	*B = (int) (0.055648*X - 0.204043*Y + 1.057311*Z); 
}


int SEEDS::RGB2LAB(const int& r, const int& g, const int& b, float* lval, float* aval, float* bval)
{
	float xVal = 0.412453 * r + 0.357580 * g + 0.180423 * b;
	float yVal = 0.212671 * r + 0.715160 * g + 0.072169 * b;
	float zVal = 0.019334 * r + 0.119193 * g + 0.950227 * b;

	xVal /= (255.0 * 0.950456);
	yVal /=  255.0;
	zVal /= (255.0 * 1.088754);

	float fY, fZ, fX;
	float lVal, aVal, bVal;
	float T = 0.008856;

	bool XT = (xVal > T);
	bool YT = (yVal > T);
	bool ZT = (zVal > T);

	fX = XT * pow(xVal,1.0/3.0) + (!XT) * (7.787 * xVal + 16.0/116.0);

	// Compute L
	float Y3 = pow(yVal,1.0/3.0); 
	fY = YT*Y3 + (!YT)*(7.787*yVal + 16.0/116.0);
	lVal  = YT * (116 * Y3 - 16.0) + (!YT)*(903.3*yVal);

	fZ = ZT*pow(zVal,1.0/3.0) + (!ZT)*(7.787*zVal + 16.0/116.0);

	// Compute a and b
	aVal = 500 * (fX - fY);
	bVal = 200 * (fY - fZ);

	*lval = lVal;
	*aval = aVal;
	*bval = bVal;

	return 1; //LAB2bin(lVal, aVal, bVal);
}



int SEEDS::RGB2HSV(const int& r, const int& g, const int& b, float* hval, float* sval, float* vval)
{
	float r_ = r / 256.0;
	float g_ = g / 256.0;
	float b_ = b / 256.0;

    float min_rgb = min(r_, min(g_, b_));
    float max_rgb = max(r_, max(g_, b_));
    float V = max_rgb;
	float H = 0.0;
	float S = 0.0;

    float delta = max_rgb - min_rgb;
	
	if ((delta > 0.0) && (max_rgb > 0.0))
	{
		S = delta / max_rgb;
	    if (max_rgb == r_)
	        H = (g_ - b_) / delta;
    	else if (max_rgb == g_)
	        H = 2 + (b_ - r_) / delta;
	    else
	        H = 4 + (r_ - g_) / delta;
	}

    H /= 6;

	/*float V = max(r_,max(g_,b_));
	float S = 0.0;
	if (V != 0.0) S = (V - min(r_, min(g_,b_)))/V;
	float H = 0.0;
	if (S != 0.0)
	{
		if (V == r_)
			H = ((g_ - b_)/6.0)/S;
		else if (V == g_)
			H = 1.0/2.0 + ((b_-r_)/6.0)/S;
		else
			H = 2.0/3.0 + ((r_-g_)/6.0)/S;
	}*/

	if (H<0.0) H += 1.0;

	*hval = H;
	*sval = S;
	*vval = V;

	if ((H<0.0) || (H>= 1.0)) printf("H: %f \n", H);

	int hbin = floor(H * nr_bins);
	int sbin = floor(S * nr_bins);
	int vbin = floor(V * nr_bins);

	//printf("%d %d %d -- %f %f %f -- bins %d %d %d\n", r, g, b, H, S, V, hbin, sbin, vbin);

	return hbin + nr_bins*sbin + nr_bins*nr_bins*vbin;
}

int SEEDS::LAB2bin(float l, float a, float b)
{
	// binning
	int bin_l = floor(l/100.1*nr_bins);

	int bin_a = floor(a/100*(nr_bins-2) + (nr_bins-2)/2 + 1);
	if (bin_a < 0) bin_a = 0;
	if (bin_a >= nr_bins) bin_a = nr_bins-1;

	int bin_b = floor(b/100*(nr_bins-2) + (nr_bins-2)/2 + 1);
	if (bin_b < 0) bin_b = 0;
	if (bin_b >= nr_bins) bin_b = nr_bins-1;

	// encoding
	return bin_l + nr_bins*bin_a + nr_bins*nr_bins*bin_b;
}


int SEEDS::RGB2LAB_special(int r, int g, int b, float* lval, float* aval, float* bval)
{
	float xVal = 0.412453 * r + 0.357580 * g + 0.180423 * b;
	float yVal = 0.212671 * r + 0.715160 * g + 0.072169 * b;
	float zVal = 0.019334 * r + 0.119193 * g + 0.950227 * b;

	xVal /= (255.0 * 0.950456);
	yVal /=  255.0;
	zVal /= (255.0 * 1.088754);

	float fY, fZ, fX;
	float lVal, aVal, bVal;
	float T = 0.008856;

	bool XT = (xVal > T);
	bool YT = (yVal > T);
	bool ZT = (zVal > T);

	fX = XT * pow(xVal,1.0/3.0) + (!XT) * (7.787 * xVal + 16/116);

	// Compute L
	float Y3 = pow(yVal,1.0/3.0); 
	fY = YT*Y3 + (!YT)*(7.787*yVal + 16/116);
	lVal  = YT * (116 * Y3 - 16.0) + (!YT)*(903.3*yVal);

	fZ = ZT*pow(zVal,1.0/3.0) + (!ZT)*(7.787*zVal + 16/116);

	// Compute a and b
	aVal = 500 * (fX - fY);
	bVal = 200 * (fY - fZ);

	*lval = lVal;
	*aval = aVal;
	*bval = bVal;

	int bin1 = 0;
	int bin2 = 0;
	int bin3 = 0;

	while (lVal > bin_cutoff1[bin1]) {
		bin1++;
	}
	while (aVal > bin_cutoff2[bin2]) {
		bin2++;
	}
	while (bVal > bin_cutoff3[bin3]) {
		bin3++;
	}

	return bin1 + nr_bins*bin2 + nr_bins*nr_bins*bin3;
}


int SEEDS::RGB2LAB_special(int r, int g, int b, int* bin_l, int* bin_a, int* bin_b)
{
	float xVal = 0.412453 * r + 0.357580 * g + 0.180423 * b;
	float yVal = 0.212671 * r + 0.715160 * g + 0.072169 * b;
	float zVal = 0.019334 * r + 0.119193 * g + 0.950227 * b;

	xVal /= (255.0 * 0.950456);
	yVal /=  255.0;
	zVal /= (255.0 * 1.088754);

	float fY, fZ, fX;
	float lVal, aVal, bVal;
	float T = 0.008856;

	bool XT = (xVal > T);
	bool YT = (yVal > T);
	bool ZT = (zVal > T);

	fX = XT * pow(xVal,1.0/3.0) + (!XT) * (7.787 * xVal + 16/116);

	// Compute L
	float Y3 = pow(yVal,1.0/3.0); 
	fY = YT*Y3 + (!YT)*(7.787*yVal + 16/116);
	lVal  = YT * (116 * Y3 - 16.0) + (!YT)*(903.3*yVal);

	fZ = ZT*pow(zVal,1.0/3.0) + (!ZT)*(7.787*zVal + 16/116);

	// Compute a and b
	aVal = 500 * (fX - fY);

	bVal = 200 * (fY - fZ);

	//*lval = lVal;
	//*aval = aVal;
	//*bval = bVal;

	int bin1 = 0;
	int bin2 = 0;
	int bin3 = 0;

	while (lVal > bin_cutoff1[bin1]) {
		bin1++;
	}
	while (aVal > bin_cutoff2[bin2]) {
		bin2++;
	}
	while (bVal > bin_cutoff3[bin3]) {
		bin3++;
	}

	*bin_l = bin1;
	*bin_a = bin2;
	*bin_b = bin3;

	return bin1 + nr_bins*bin2 + nr_bins*nr_bins*bin3;
}

int SEEDS::threebythree(int x, int y, int label)
{
	int count = 0;
	if (labels[seeds_top_level][(y-1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x+1]==label) count++;

	if (labels[seeds_top_level][y*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][y*img_width+x+1]==label) count++;

	if (labels[seeds_top_level][(y+1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x+1]==label) count++;

	return count;
}

int SEEDS::threebyfour(int x, int y, int label)
{
	int count = 0;
	if (labels[seeds_top_level][(y-1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x+1]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x+2]==label) count++;

	if (labels[seeds_top_level][y*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][y*img_width+x+2]==label) count++;

	if (labels[seeds_top_level][(y+1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x+1]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x+2]==label) count++;

	return count;
}

int SEEDS::fourbythree(int x, int y, int label)
{
	int count = 0;
	if (labels[seeds_top_level][(y-1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y-1)*img_width+x+1]==label) count++;

	if (labels[seeds_top_level][y*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][y*img_width+x+2]==label) count++;

	if (labels[seeds_top_level][(y+1)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y+1)*img_width+x+2]==label) count++;

	if (labels[seeds_top_level][(y+2)*img_width+x-1]==label) count++;
	if (labels[seeds_top_level][(y+2)*img_width+x]==label) count++;
	if (labels[seeds_top_level][(y+2)*img_width+x+1]==label) count++;

	return count;
}


float SEEDS::geometric_distance(int label1, int label2)
{
	float dx = ((float)x_position[label1]/T[seeds_top_level][label1] - (float)x_position[label2]/T[seeds_top_level][label2]);
	float dy = ((float)y_position[label1]/T[seeds_top_level][label1] - (float)y_position[label2]/T[seeds_top_level][label2]);
	return sqrt(dx*dx + dy*dy);
}



float SEEDS::intersection(int level1, int label1, int level2, int label2)
{
    float intersect = 0.0;
	
	for (int n=0; n<histogram_size; n++)
	{
		intersect += min((float)histogram[level1][label1][n]/T[level1][label1], (float)histogram[level2][label2][n]/T[level2][label2]);
	}

	return intersect;
}





bool SEEDS::check_split(int a11, int a12, int a13, int a21, int a22, int a23, int a31, int a32, int a33, bool horizontal, bool forward)
{
	if (horizontal)
	{
		if (forward)
		{
			if ((a22 != a21) && (a22 == a12) && (a22 == a32)) return true;
			if ((a22 != a11) && (a22 == a12) && (a22 == a21)) return true; 
			if ((a22 != a31) && (a22 == a32) && (a22 == a21)) return true;
		}
		else /*backward*/
		{
			if ((a22 != a23) && (a22 == a12) && (a22 == a32)) return true;
			if ((a22 != a13) && (a22 == a12) && (a22 == a23)) return true; 
			if ((a22 != a33) && (a22 == a32) && (a22 == a23)) return true;
		}
	}
	else /*vertical*/
	{
		if (forward)
		{
			if ((a22 != a12) && (a22 == a21) && (a22 == a23)) return true;
			if ((a22 != a11) && (a22 == a21) && (a22 == a12)) return true; 
			if ((a22 != a13) && (a22 == a23) && (a22 == a12)) return true;
		}
		else /*backward*/
		{
			if ((a22 != a32) && (a22 == a21) && (a22 == a23)) return true;
			if ((a22 != a31) && (a22 == a21) && (a22 == a32)) return true; 
			if ((a22 != a33) && (a22 == a23) && (a22 == a32)) return true;
		}
	}

	return false;
}


int SEEDS::count_superpixels()
{
	int* count_labels = new int[nr_labels[seeds_top_level]];
	for (int i=0; i<nr_labels[seeds_top_level]; i++)
		count_labels[i] = 0;
	for (int i=0; i<img_width*img_height; i++)
		count_labels[labels[seeds_top_level][i]] = 1;
	int count = 0;
	for (int i=0; i<nr_labels[seeds_top_level]; i++)
		count += count_labels[i];

	delete count_labels;
	return count;
}


void SEEDS::SaveLabels_Text(string filename) 
{
	int sz = img_width*img_height;

	ofstream outfile;
	outfile.open(filename.c_str());
        int i = 0;
        for( int h = 0; h < img_height; h++ )
          {
            for( int w = 0; w < img_width -1; w++ )
              {
                outfile << labels[seeds_top_level][i] << " ";
                i++;
              }
            outfile << labels[seeds_top_level][i] << endl;
            i++;
          }
	outfile.close();
}


void SEEDS::compute_mean_map()
{
	means = new UINT[img_width*img_height];

	for (int i=0; i<img_width*img_height; i++)
	{
		int label = labels[seeds_top_level][i];
		float L = 100.0 * ((float) L_channel[label]) / T[seeds_top_level][label];
		float a = 255.0 * ((float) A_channel[label]) / T[seeds_top_level][label] - 128.0;
		float b = 255.0 * ((float) B_channel[label]) / T[seeds_top_level][label] - 128.0;
		int R, G, B;
		LAB2RGB(L, a, b, &R, &G, &B);
		means[i] = B | (G << 8) | (R << 16);
	}
}












